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In van der Waals bonded or rotationally disordered mnltilayer stacks of two- 
dimensional (2D) materials, the electronic states remain tightly confined within 
individual 2D layers. As a result, electron-phonon interactions occur primar¬ 
ily within layers and interlayer electrical conductivities are low. In addition, 
strong covalent in-plane intralayer bonding combined with weak van der Waals 
interlayer bonding results in weak phonon-mediated thermal coupling between 
the layers. We demonstrate here, however, that Coulomb interactions between 
electrons in different layers of multilayer epitaxial graphene provide an impor¬ 
tant mechanism for interlayer thermal transport even though all electronic states 
are strongly confined within individual 2D layers. This effect is manifested in 
the relaxation dynamics of hot carriers in ultrafast time-resolved terahertz spec¬ 
troscopy. We develop a theory of interlayer Coulomb coupling containing no free 
parameters that accounts for the experimentally observed trends in hot-carrier 
dynamics as temperature and the number of layers is varied. 
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Introduction 


The dynamics of electrons in atomic-layer 2D electron systems like graphene is a snbject of 
considerable cnrrent interest, partly becanse of its relevance to a wide variety of potential 
electronic and optoelectronic device applications. Many proposed and prototype devices em¬ 
ploy stacks composed of many 2D electronic material layers. Examples of mnltilayer systems 
inclnde mnltilayer epitaxial graphene (MEG),^ van der Waals bonded layered sheets,® 
transition metal dichalcogenides^^ and others.® In these strnctnres the interactions between 
electrons in different layers in the stack becomes a snbject of key importance. One important 
property is that phonon-mediated interlayer thermal conpling is weak relative to that in bnlk 
3D materials. In the example of MEG, rotational stacking arrangements deconple electronic 
states localized in different 2D layers.®®' As a resnlt, phonon-mediated interlayer thermal 
conpling in MEG is strongly rednced relative to typical bnlk behavior in 3D materials. 

The qnestion thns arises as to what other mechanisms can contribnte to thermal eqnili- 
bration between different layers. We consider this qnestion here in the context of hot carrier 
dynamics. If electrons are heated in one layer (e.g. by optical excitation or electrical injec¬ 
tion), they will normally cool to the lattice temperatnre by optical phonon emission at high 
carrier energies, and by aconstic phonon emission at low carrier energies. For graphene, it 
is well established that electron cooling by aconstic phonons is very efficient in highly doped 
layers.'®®®' The sitnation is qnite different, however, for lightly doped or nearly nentral 
graphene in which a small joint-density-of-states for electronic transitions combines with 
a small acoustic phonon energy at typical scattering wavevectors to diminish the acoustic 
phonon cooling power.'®®' (The cooling power in this limit can, however, be substantially 
enhanced by disorder, because it relaxes momentum conservation limits'®®®' on allowed 
processes.) 

An interesting case thus arises when a multilayer stack contains both highly doped (HD) 
and nearly neutral lightly doped (LD) graphene layers. This is exactly the situation that 
occurs in MEG grown on the G-face of SiG,'®' and is likely to be relevant to gated multilayer 


3 








2D systems due to interlayer screeningJ ^^ * ^^! If electrons are heated in the multilayer struc¬ 
ture, then acoustic-phonon-mediated cooling would result in the rapid buildup of a thermal 
gradient between the HD and LD layers; the HD layers would quickly approach the lattice 
temperature, while carriers in the LD layers would remain hot. Eventually, of course, ther¬ 
mal equilibrium would be restored, as thermal energy flows from the LD to the HD layers. 
The HD layers can be a heat sink for the LD layers if there is an effective interlayer energy 
transfer mechanism. 

In the following, we show that Coulomb scattering between electrons in LD and HD 
layers of MEG can provide an efficient means for interlayer thermal coupling, and provide 
an alternate mechanism for cooling of hot electrons in the LD layers that acts in parallel 
with acoustic-phonon-mediated intralayer cooling. This process is illustrated schematically 
in Figure [^. We note that interlayer thermal coupling via Coulomb scattering has been 
considered recently in the context of 2D electron gases in transport devices.^ We begin by 
outlining a heuristic analytical model calculation for a pair of graphene layers to establish 
the magnitude of the effect relative to acoustic-phonon-mediated intralayer cooling. This 
simple calculation establishes the signihcance of the effect; we then discuss how hot-carrier 
cooling in multilayer systems is accessible experimentally by ultrafast time-resolved terahertz 
(THz) spectroscopy and ultrafast infrared (IR) pump-probe spectroscopy measurements. 
Following a discussion of the features of the data that point to interlayer energy transfer, we 
present the details of a theory of interlayer energy transfer via screened Coulomb interactions. 
The calculated cooling powers imply asymptotic cooling times on the sub-nanosecond scale. 
We show that the calculated dynamics and trends with lattice temperature and number of 
epitaxial graphene layers are fully consistent with the experimental results, without the need 
for any htting parameters. 
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Results 


Interlayer Coulombic energy transfer heuristics 

We first ask whether or not interlayer Coulomb coupling can potentially dominate over 
acoustic phonon coolin^^ and disorder-assisted electron-phonon (supercollision) coolin^^ 
in multilayer graphene samples. A simple comparison of the cooling powers of the different 
mechanisms suggests that the answer is yes. The low-temperature cooling power of 
interlayer Coulombic energy transfer between a hot LD and a cold HD graphene layer is (see 
Supplementary Notes 2-3): 


V ^bT J 


( 1 ) 


where Ui = 2E-F^i/(irh^Vp) is the density of states at the Fermi level. Equation is plotted 
in Figure]^ as a function of electron temperature for Fermi level Ef,hd = 300 meV in the 
HD layer and various Fermi levels Ef,ld in the LD layer; even at very low carrier density 
in the LD layer, the cooling power is quite substantial. We can compare Equation with 
the cooling powers of both acoustic phonon cooling Q®' oc T and disorder-assisted electron- 
phonon (supercollision) cooling Q®® oc of an isolated LD layer. The ratio of cooling 
powers is plotted in Figure]^ as a function of electron temperature for Ef,hd = 300 

meV and various values of Ef,ld- As it can be expected, acoustic phonon cooling is very 
inefficient in graphene with very low carrier density. The ratio of cooling powers Q®® is 
also plotted in Figure as a function of electron temperature (above the Bloch-Griineisen 
temperature Tbg ~ 5 K) for Ef,hd = 300 meV, Ef,ld = 10 meV (typical for C-face MEG 
on SiG) and various values of the low-density-layer disorder mean free path. It is apparent 
that for high quality graphene such as G-face MEG on SiG,^ interlayer Goulombic energy 
transfer can dominate for a wide range of electron temperatures and sample characteristics. 
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Multilayer epitaxial graphene 


In the main body of this paper, we investigate the physics of interlayer Conlomb coupling 
in MEG, grown on the C-face of single-crystal 4H-SiC(0001) substrates by thermal decom¬ 
position of Si atoms.l^ An important feature of this material is that it largely preserves 
distinct graphene-like electronic properties because of unique rotational stacking, which sup¬ 
presses hybridization between low energy electronic states localized in neighboring planes of 
carbon atoms.^^^ MEG is doped by electron transfer from the interface with the support¬ 
ing SiG substrate and the induced n-type carrier-density prohle falls off rapidly with layer 
moving away from the substrate (see the inset of Figure]^). We will refer to the few layers 
close to the SiG substrate, which have large carrier densities of nno ^ lO^^cm”^ as deter¬ 
mined from high-resolution angle-resolved photoemission spectroscopy (ARPES), scanning 
tunneling spectroscopy (STS), electronic transport, and ultrafast optical spectroscopy mea- 
surements,^^^2III126l2ll high-density (HD) layers, and to those further away, whose carrier 
densities drop quickly*^ to the range of uld ^ 10^°cm“^ as determined from STS, electronic 
transport, and magneto-optical spectroscopy measurements,II^ as low-density (LD) layers. 
The formation of local spatial charge inhomogeneities due to small amounts of disorder, 
impurities, or surface corrugation of the SiG substrate could explain the non-zero carrier 
density measured in the top layers of MEG.^^^^ 

Electronic cooling in multilayer epitaxial graphene 

When a MEG sample is illuminated with a short optical pulse, electrons are excited to 
high energies, leaving behind unoccupied states or holes. Due to strong intralayer carrier- 
carrier scattering, these hot carriers thermalize with the background of cold carriers within 
~ 50 fs,l3H2Hl forming two separate non-equilibrium Fermi-Dirac distributions for electrons 
and for holes. The electron and the hole quasi-Fermi levels subsequently merge within 
~ 100 — 200 fs,l333El establishing a uniform electron liquid within each layer i characterized 
by an elevated electron temperature, T,. As the electron liquid cools, each layer’s electron 
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temperature approaches the equilibrium lattice temperature, Tl. It is generally accepted 
that initial fast cooling occurs in the hrst few picoseconds via the emission of energetic 
optical phonons, and that this process becomes increasingly inefficient as the electron energy 
falls below the relatively high optical phonon energy {hoOop ~ 200 meV^^. In the final stage 
of relaxation, low energy electronic cooling in an isolated layer would proceed by the much 
slower emission of acoustic phonons. Recent findings suggest that that the low energy phonon 
spectra of multilayer graphene systems are sensitive to the pattern of relative orientations.*^ 
Although this property will transfer a corresponding sensitivity to phonon cooling powers, 
relative orientations will not influence the interlayer Coulombic energy transfer mechanism 
explored here. Prior theoretical work has found that the rate of acoustic phonon cooling 
in disorder-free single-layer graphene is very strongly dependent on the carrier density, with 
cooling times ranging from tens of picoseconds for doping densities of ~ lO^^cm”^ to tens of 
nanoseconds for doping densities of ~ The strong carrier density dependence 

of the acoustic phonon emission implies that the LD and HD layers of MEG will exhibit 
very differing cooling rates following ultrafast optical excitation, leading to the buildup of a 
thermal gradient, which triggers an interlayer energy transfer. 

We have applied two different experimental techniques to probe the dynamics of the 
interlayer energy transfer in MEG. The first is ultrafast time-resolved THz spectroscopy, 
which is a powerful tool for investigating the real time relaxation dynamics of photoexcited 
carriers, because it is sensitive to both the number of carriers and their distribution in 
energyDue to the large number of layers in our MEG samples and the rapid decrease 
of carrier density with layer number, the measured differential THz transmission signal is 
dominated by the dynamic THz response of the many LD layers and reveals the cooling of 
the LD layers due to their coupling to the HD layers. The second is ultrafast degenerate IR 
pump-probe spectroscopy, in which we optically inject hot carriers selectively into the LD 
layers and then observe directly the transfer of heat to the electron liquid in the most highly 
doped HD layer.*^^*^ 
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Ultrafast time-resolved THz spectroscopy 

We first consider the ultrafast time-resolved THz spectroscopy experiments. We report 
measurements on a series of MEG samples ranging from ~ 3 to ~ 63 layers in an ultrafast 
optical-pump THz-probe set-upA^*^ as illustrated schematically in Figure [^. Our laser 
system consists of a Ti:Sapphire oscillator (Mira 900-F, Coherent) followed by a Ti;Sapphire 
regenerative ampliher (RegA 9050, Coherent) and produces ultrafast optical pulses with a 
center wavelength of 800 nm, a pulse width of ~ 60 fs and a repetition rate of 250 kHz. 
A portion of the laser beam is quasi-collimated at the sample position with an intensity 
spot size diameter of ~ 1600 /rm, and optically injects hot carriers in the MEG samples. 
A second portion of the laser beam is used to generate a single-cycle THz pulse in a low 
temperature grown GaAs photoconductive emitter (Tera-SED 3/4, Gigaoptics)^^^^ and the 
emitted broadband THz radiation is focused on the MEG sample with an intensity spot 
size diameter of ~ 500 pm to probe the dynamic THz response. The transmitted portion 
of the THz probe is detected by using time-domain electro-optic sampling in a 1 mm thick 
ZnTe crystall^^® and a pair of balanced Si photodiodes. The electrical signal is modulated 
by a mechanical chopper, placed in either the optical pump or the THz probe arm, and 
recorded by using a conventional lock-in ampliher data acquisition technique. The MEG 
sample is mounted inside a liquid helium continuous how cryostat (ST-100, Janis), so that 
the substrate temperature can be varied from 10 K to 300 K. The time delays between the 
optical pump, the THz probe and the sampling pulse are controlled by two motorized stages. 
All THz optics is surrounded by an enclosure purged with purihed nitrogen gas to minimize 
water vapor absorption. The detection bandwidth of the system is in the range of ~ 0.2 — 2.5 
THz and the temporal resolution of the measurements is limited by the duration of the THz 
probe pulse to the sub-picosecond timescale. The experimental error is due primarily to long¬ 
term drift of the optomechanical components and the ultrafast Ti:Sapphire laser system, and 
is estimated not to exceed ~ 5%. 

As a hrst experimental approach, we measure the diherential change in the THz probe 






pulse transmission through the MEG sample due to photoexcitation. The THz probe held 
is recorded in the time-domain, and it is later numerically Fourier transformed to obtain the 
frequency spectrum. Figure shows the differential THz transmission spectra normalized 
to the THz transmission without photoexcitation, At{u)/t{u), for a few different THz probe 
delays after the optical pump for a MEG sample with ~ 63 layers. From the Tinkham 
formula for the transmission through a thin conducting him on a transparent substrate,*^ 
the At{uj)/t{u) signal can be directly related to the photoinduced change in the complex 
sheet conductivity of MEG, A(t(c(;), through the expression: 

. , . nsuh + ^ivac At{u) 

Aa{uj) ^ - X , , , (2) 

^ ^ ho t{u) ^ ^ 

where Ugub and Uvac are the THz refractive indices of the SiG substrate and the environ¬ 
ment, respectively, and tjq is the impedance of free space. The THz conductivity of MEG 
is a summation of the THz conductivities of the individual epitaxial graphene layers in the 
MEG stack, because the layers are electronically decoupled.*^*^ Additional THz transmission 
spectra, similar to the ones in Figure [^, but for variable substrate temperature and variable 
pump huence, are shown in Supplementary Fig. 1-2 and Supplementary Note 1. We note 
that the normalized diherential THz transmission spectra are remarkably dispersionless in 
the detectable frequency range under all experimental conditions; this justihes the applica¬ 
tion of a simpler data acquisition scheme in which we record the normalized differential THz 
transmission only at the peak of the THz probe pulse. 

As a second experimental approach, we keep the delay of the sampling pulse fixed at the 
peak of the THz probe pulse, and we scan the pump-probe delay to map out the relaxation 
dynamics of the photoexcited carriers. Since the carrier-carrier scattering time in graphene 
is much shorter than the temporal duration of the THz probe, we study the relaxation of 
the THz transmission (or the THz conductivity) change induced by the optical pump in the 
limit, where it is determined by collective electronic cooling dynamics. Figure [^-d show 


9 




the normalized differential THz transmission at the peak of the THz probe pulse, At/t, 
as a function of pump-probe delay for variable substrate temperature for the same MEG 
sample with ~ 63 layers. At time zero, the optical pump photoexcites carriers in the MEG 
sample resulting in an overall increase of the THz conductivity and hence THz absorption, 
as manifested in a negative differential THz transmission. In graphene with very low doping, 
the increase of the electron temperature leads primarily to larger electron occupation in the 
conduction band and a corresponding net increase of the THz conductivity, consistent with 
our interpretation that the measured dynamic THz response is dominated by the hot carriers 
in the many LD layers of MEG.^^^^The differential THz transmission reaches its maximum 
magnitude within ~ 1 ps, with the rise time being limited mainly by the temporal duration of 
the THz probe. The differential THz transmission subsequently recovers as the thermalized 
hot carriers cool to the substrate temperature with relaxation times ranging from a few 
picoseconds at room temperature to hundreds of picoseconds at cryogenic temperatures. 
The secondary decrease in the differential THz transmission at ~ 7 ps is due to a round-trip 
reflection of the optical pump inside the substrate that photoexcites additional carriers. 

We perform phenomenological fits to the experimental data in Figure [^-d, and we dis¬ 
cover that the differential THz transmission evolves from a faster mono-exponential decay 
at room temperature to a slower bi-exponential decay at cryogenic temperatures. As we 
explain below, the slow electronic cooling at low substrate temperatures is controlled by 
interlayer Goulombic energy transfer between the LD and HD layers. A summary of the 
extracted carrier relaxation times as a function of substrate temperature for a few different 
pump ffuences is presented in Figure [^. We observe that the relaxation times are largely in¬ 
dependent of the pump ffuence (or the initial electron temperature) except at high substrate 
temperatures. The slight increase in the relaxation times at the highest pump ffuence can be 
attributed to heating of the HD layers above the substrate temperature, which decreases the 
rate of the energy transfer from the LD layers. Similarly, the energy transfer between layers 
becomes less efficient at high substrate temperatures, at which the difference between the 
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electron temperatures in different layers is small. As a consequence, the contribution of the 
interlayer Coulomb coupling to electronic cooling is diminished at substrate temperatures 
above ~ 200 K as evidenced from the hts. 

We repeat identical experiments and analysis for a second MEG sample with ~ 35 layers 
and the corresponding summary of the extracted carrier relaxation times for the best hts 
are presented in Figure]^. Qualitatively, the THz carrier dynamics for the 35-layer sample 
mirror those for the 63-layer sample by exhibiting a transition from a faster mono-exponential 
to a slower bi-exponential decay as the substrate temperature is decreased. Similarly, the 
relaxation times are independent of the pump huence except at high substrate temperatures. 
Further inspection shows that the long relaxation times become up to a few times shorter, 
when the number of epitaxial graphene layers is nearly halved, which indicates the presence of 
interlayer interaction. We show below that because of the range dependence of the Coulomb 
scattering processes, the addition of more LD layers slows their collective electronic cooling 
via coupling to the HD layers. 

To underscore the profound influence of interlayer energy transfer on the electronic cooling 
in MEG, we next study the limiting case of MEG with all HD layers. We again perform 
identical experiments and analysis for a third MEG sample with only ~ 3 layers. Figure 
shows the normalized differential THz transmission at the peak of the THz probe pulse, 
At/t, as a function of pump-probe delay for variable substrate temperature. First, we note 
that the differential THz transmission for MEG with all HD layers is positive, which has 
been previously phenomenologically attributed to enhanced carrier scattering as the electron 
temperature is elevated.®^ Second, we observe that the THz carrier dynamics are much 
faster and completely independent of the substrate temperature, because there is practically 
very little or no interlayer energy transfer. The differential THz transmission is £t very 
well by a phenomenological mono-exponential decay at all temperatures (again accounting 
for the substrate reflection of the optical pump) and the summary of the extracted carrier 
relaxation times as a function of substrate temperature for a few different pump fluences 
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is presented in Figure it These relaxation dynamics are inconsistent with the disorder- 
assisted electron-phonon (supercollision) cooling mechanism (see Supplementary Fig. 5-10 
and Supplementary Note 6), but they could be attributed to the optical and acoustic phonon 
cooling mechanisms of hot carriers. 

Ultrafast degenerate IR pump-probe spectroscopy 

To further support the existence of interlayer energy transfer in MEG, we devise another 
experiment using ultrafast degenerate IR pump-probe spectroscopy^^SEll in which we selec¬ 
tively photoexcite hot electrons in all layers of MEG except the hrst HD layer nearest to the 
SiG substrate. We then observe the cooling of these hot electrons via interlayer Goulomb 
coupling to the cold electrons in the hrst HD layer. Our laser system consists again of a 
Ti:Sapphire oscillator and ampliher followed by an optical parametric ampliher (OPA 9850, 
Goherent) and produces ultrafast optical pulses with a center wavelength tuned to 1.8 pm. 
The OPA beam is hltered through a 10 nm bandpass hlter centered at 1.8 pm and is split 
into a pump and a probe beams that are both focused on the MEG sample with an intensity 
spot size diameter of ~ 50 pm. The transmitted portion of the probe is detected by using 
a grating spectrometer and an InGaAs photodetector in conjunction with a conventional 
lock-in ampliher data acquisition technique. The MEG sample is mounted inside a cryostat 
(ST-100, Janis) and the substrate temperature is held at 10 K. The experimental error is 
again estimated not to exceed ~ 5%. 

In this experimental approach, both the pump and the probe photon energies (hca ~ 690 
meV) are chosen to be slightly smaller than twice the Fermi level of the hrst HD layer 
{Ep ^ 360 meV), but larger than twice the Fermi levels of all other layers in MEG.^^^^ 
Thus, the pump selectively injects hot electrons in all layers of MEG except the hrst HD 
layer, in which interband absorption is Pauli blocked as illustrated schematically in the 
inset of Figure The probe diherential transmission due to photoexcitation has a positive 
contribution from the hot electrons in all layers of MEG except the hrst HD layer. The hrst 
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HD layer gives no contribution, when the carriers in that layer remain unexcited. However, 
interlayer Coulombic energy transfer can heat this layer, which results in a distinct negative 
contribution to the differential transmission. 

Figure 1^ shows the probe differential transmission normalized to the probe transmission 
without photoexcitation, AT/T, as a function of pump-probe delay for the MEG sample with 
~ 63 layers. We observe that immediately after photoexcitation the differential transmission 
is positive, arising from the hot electrons that are directly injected in the top layers. Shortly 
after that, the differential transmission becomes negative and reaches its minimum value 
within ~ 1 ps. This sign change demonstrates the existence of an efficient interlayer energy 
transfer, in which the cold electrons in the hrst HD layer act as a heat sink for the hot 
electrons in the top layers. The rapidly rising electron temperature in the hrst HD layer has 
a dominant negative contribution to the differential transmission. As delay time increases, 
the electrons in the hrst HD layer cool much faster than the electrons in the top layers due to 
the sharply increasing rate of acoustic phonon emission with carrier density. This results in a 
second sign change in the diherential transmission at ~ 20 ps, when the electron temperature 
in the hrst HD layer approaches the equilibrium lattice temperature. By comparing to the 
THz carrier dynamics in Figure [^, we note that it takes slightly longer for the HD layer 
to cool due to the extra heat from the top LD layers. At that point, the optical and 
acoustic phonon cooling rate in the hrst HD layer balances the interlayer energy transfer 
rate from the top layers. Electronic cooling in the top layers of MEG via the interlayer 
Goulomb coupling survives on a timescale exceeding a hundred picoseconds (limited by the 
experimental signal-to-noise ratio), which is consistent with that observed in the THz carrier 
dynamics in Figure |^-d. Now, we turn our attention to a detailed presentation of the 
theory of hot-carrier equilibration based on interlayer energy transfer via screened Goulomb 
interactions. 
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Interlayer Coulombic energy transfer theory 

Non-equilibrium electrons in graphene have been shown to thermalize within a layer on an 
ultrafast timescale on the order of tens of femtoseconds.®® We assume that this property 
holds also in MEG, and that it leads to pseudo-equilibrium electronic states with well defined 
temperatures Tj in layer i, but we allow for the possibility of differences in temperature 
between layers that survive to longer timescales. The multilayer temperature dynamics are 
described by a set of coupled non-linear first order differential equations: 

s,r, = (Ti) + eS(r.. I)) j /Ci, (3) 

where C* = drSi is the heat capacity and Si the energy density of electrons in the i’th layer. 
The rate of change of energy density in layer i is determined by the sum of two processes: 
energy loss to the lattice via electron-phonon scattering at rate and energy transfer 

from other layers j to i via interlayer electron-electron scattering at rate Qfy Both of these 
mechanisms depend strongly on the carrier density. The four most highly doped layers in 
the 63-layer MEG sample have Fermi energies measured to be 360, 218, 140, and 93 meV, 
respectively.® A simple Thomas-Fermi model® is able to account semi-quantitatively for 
the monotonic decrease in carrier density with separation from the substrate in multilayer 
graphene systems. Electronic transport and magneto-optical spectroscopy measurements^^® 
of top LD layers suggest local carrier density fluctuations in these LD layers that satisfy 

^LD < 10^°cm“T 

Recent ultrafast optical spectroscopy experiments® have found that the hot carriers in 
the HD layers of MEG quickly relax to the lattice temperature with equilibration times on 
the order of a few picoseconds, in good agreement with the experiments reported here. The¬ 
oretical calculations neglecting interlayer thermal coupling (Qd —)■ 0) have obtained order of 
magnitude agreement with these measurements.® On the other hand, the same approxima¬ 
tion applied to the LD layers erroneously predicts thermal equilibration times on the order of 
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several nanoseconds, in sharp disagreement with the experiments reported here. This was the 
initial impetus for our examination of the interlayer Coulombic energy transfer mechanism. 
According to Ref. HD the acoustic phonon cooling power is proportional to the square of the 
carrier density, n^. Allowing for disorder-assisted electron-phonon (supercollision) scattering 
changes this dependence to In the HD layers, optical and acoustic phonon emission is 
likely to provide the dominant cooling pathway for hot carriers.!^ In the LD layers, however, 
we suggest that it plays a more subsidiary role by keeping the electron temperature in the HD 
layers pinned to the lattice temperature (Thd = Rl), while they act as a Coulomb-coupled 
heat sink for the LD layers (see Supplementary Fig. 4 and Supplementary Note 5). 

Although the true electron density prohle across the many layers of MEG is expected to 
decrease smoothly, it is convenient to make a sharp distinction between highly doped (HD) 
and lightly doped (LD) layers. We predict that the asymptotic temperature dynamics of 
LD layers in MEG are effectively governed by Equation with 0- theoretical 

analysis, we will denote the four most highly doped layers near the substrate, as HD layers. 
Quantitatively, this cutoff is suggested from the disorder-free theory of acoustic phonon 
cooling,^ which in combination with Equation 4 in Supplementary Note 2, implies that 
the hot-electron distribution in layers i > 4 transfers energy via the interlayer Coulomb 
interaction to the j HD layers {j < i) faster than it loses energy to the lattice via acoustic 
phonon emission. More precisely, we hnd that Qfj) ~ temperatures T* > 50 

K. In all calculations described below, we use the values Ef,i -4 = 360, 218, 140, and 93 meV 
for the Fermi levels of the hrst four HD layers. These have been measured explicitly for the 
63-layer MEG sample and are expected to be good estimates for the 35-layer MEG sample. 

The rate of Coulombic energy transfer between two layers is given by a Fermi golden-rule 
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expression: 


2« = 

ki,k! kj.k'. 

■’1 J’ j 

X (^k^z+k /,ki+k, <^(£^k / + £^k^z ^k, ^'ki 


X (^/kj (l - /k') /ki (1 - /k:) - /kj' (1 - /kj) /ki' (1 - fkS) , 


( 4 ) 


where (i, j) are layer indices, and k is a collective index which implies, in addition to wavevec- 
tor, the spin, the valley, and the band index labels required to specify single electron states 
in graphene’s low energy Dirac model.HSI As mentioned above, intra-layer electron thermal- 
ization®E3 jg much faster than inter-layer energy transfer, and this allows us to describe 
electronic state occupations with a quasi-equilibrium Fermi distribution /k = f{ek,^,T). In 
the random phase approximation (RPA): 




( 5 ) 


where q = |kj/ — kj|, e = — Skj, and (g, e) is the screened electron-electron interaction 

between layers i and j (see below). When screening is neglected, e) —)■ Vij{q) = 

2716“^ exp{—qdij)/Kq is the 2D Fourier transform of the bare Coulomb interaction between 
two electrons with interlayer separation dij and k = 5.5 to account for the presence of MEG 
at the surface of the SiC substrate. The factors in parenthesis are the well-known form 
factors that account for the sub-lattice spinor dependence of graphene vr-band plane-wave 
matrix elements. The indices a and fd are equal to 1 and —1 for conduction and valence band 
states respectively. The Kronecker delta’s in Equation explicitly exhibit the property that 
continuum model interactions are independent of spin (a) and valley (r). Using Equation]^ 
and comparing with the expression for graphene’s non-interacting density-density response 
function, xidi^i T), we hnally obtain the following compact expression which is suitable for 
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numerical evaluation: 


fe roo _^ 

Qg=-/ 

J-oo - 

x[nB{f\i^/Ti) - n^{hhj/Tj)] 

where n-B^x) = l/(exp(x) — 1) and /cb = 1 throughout. We use this expression below to 
calculate interlayer energy transfer rates. A central quantity in the theoretical formulation 
of the many-body effects of Dirac fermions is the dynamical polarizability tensor Xj(q', ca, Tj) 
for the z’th layer at temperature Tj. This is dehned through the one-body non-interacting 
Green’s functions.^S The density-density response function of the doped 2D Dirac electron 
model was hrst considered by Shung^SI at zero-temperature as a step toward a theory of 
collective excitations in graphite. The Dirac electron expression a;,Tj) at finite tem¬ 
perature has been recently considered.®® Before proceeding with the rate calculations, 
however, we must hrst discuss the approximation we use for the screened interlayer Coulomb 
interaction. 

In the RPA, the bare interlayer Coulomb interaction is screened® by the potential pro¬ 
duced by self-consistently readjusted charge density. For a general multilayer system this 
implies that:® 

(7) 

where v and are matrices which describe bare and screened potentials in one layer [i) 
to an external charge in another layer (j). Screening complicates energy cooling dynamics, 
because it causes the energy transfer rate between a particular pair of layers to depend, 
through y;(g,ci;,T), on the temperatures in all layers. However, the low carrier densities 
in the layers far from the substrate^ motivate a simplifying approximation in which their 
contributions to screening are neglected. By comparing this approximation with the full 
expression, we hnd that this simplihcation is justihed in the regions of phase space {q, u) 
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important for low temperature energy transfer. 

The relative ability of intraband excitations in LD and HD layers to screen the interaction 
can be established by examining the ratio of their density response functions in the static 
limit: Xhd(?)/xld(?) = /’^ld • In our MEG samples this ratio varies within ~ 8 — 31 in 

the important regions, and suggests that the leading order charge polarization responsible for 
screening the MEG interlayer Goulomb interaction can be approximated without including 
the contribution from the LD layers. We note that the Bose factors in Equation limit 
transfer energies to hu < k-^T. We also note that the factor exp(—gdjj) in the interlayer 
Goulomb interaction limits the important wavevector transfers to q < 1/dij. The conditions 
(setting h and /cb —t 1), 

q< 1/dij, (8) 

apply equally well to Goulomb-mediated interlayer energy and momentum transfer in any 
type of multilayer 2D electron system. 

Setting Xjidy^yTj) fo zero for all but the four HD layers nearest the substrate, and 
letting the separation distance between the four HD layers also go to zero (relative to the 
generally much larger spacing between HD and LD layers in MEG), we obtain the following 
approximate expression: 

v^^{q,u)ij Vij{q)/e^^^{q,u), (9) 

where 

e^’^^(g,a;) = 1 - ^ Xj(g,u;,Tj). (10) 

iSHD 

In the calculations reported on below j was summed over the four HD layers. 

Since cooling powers between LD layers are typically 1 to 2 orders of magnitude larger 
than between LD and HD layers^ (Qld ld' ^ 2ld hd); energy transfer between pairs of LD 
layers cannot be ignored in the calculations. Poles in the screened interaction described by 
Equation at plasmon frequencies greatly enhance the interlayer quasi-particle scattering 
rate. This effect does not contribute to energy transfer between HD and LD layers, because 
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the plasmon modes then lie at higher frequencies than can be excited in low temperature 
HD layers; the plasmon poles reside above the mtra-band particle-hole continuum of HD 
layers, i.e. ti;pi(g) > v-pq. However, at temperatures comparable to the Fermi temperature, 
mfer-band particle-hole excitations are no longer Pauli blocked at any frequency, and these 
excitations can take advantage of the plasmon poles in the screened interaction.^^^^ It is 
the small Fermi temperatures (Tp oc ^/n) of the LD layers that dramatically increases the 
Coulomb-coupling amongst the group of LD layers. 

As LD layers closer to the substrate cool, they absorb energy from and cool more distant 
LD layers. This effect results in a collective cooling state in which the electron temperatures 
of all LD layers relax to the equilibrium lattice temperature nearly uniformly, and motivates 
an approximation with a single collective LD layer temperature, Tc{t). In employing this 
approximation our goal is to establish that Coulomb scattering is a relevant energy transfer 
process up to quite high temperatures. A more detailed calculation in which the temper¬ 
ature of each LD layer is allowed to vary independently would be warranted if the charge 
density prohle in the multilayer system was accurately known. Such a calculation might in 
any event not achieve greater accuracy since the calculations of electron-electron scattering 
amplitudes can only be performed approximately. The collective cooling state model we em¬ 
ploy has the advantage that it requires fewer computationally burdensome hnite-temperature 
dynamic polarizability calculations. The collective cooling dynamics model of the collective 
temperature is given by: 

QfjiTuD = Tl, Tld ^ Te, d,,) j /IVCld, (11) 

VieLD isHD / 

where N is the total number of LD layers. We have summed over the energy transfer rate 
between all LD-HD pairs of layers, appealing to strong electron-phonon coupling to keep 
the HD layers temperature at Tp and to strong ld' keep the LD layers at a com¬ 
mon temperature Tc{t). Figure |^-b compares the calculated collective thermal relaxation 
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dynamics Tc{t) in 30-layer and 60-layer MEG samples for lattice temperatnre of 10 K and 
50 K. For both lattice temperatnres, thicker MEG samples cool more slowly. This trend can 
be nnderstood by noting that each additional LD layer contribntes an eqnal amonnt to the 
combined LD layers heat capacity, whereas the energy transfer rate to the HD layers falls off 
npon moving fnrther away from the snbstrate (see Snpplementary Fig. 3 and Snpplementary 
Note 4). 

If we linearize the collective temperatnre of the LD layers abont the eqnilibrinm lattice 
temperatnre, 6T = Tc{t) — Ti^, we hnd that the interlayer energy transfer rate can be written 
as: 



4Tl sinh^(ha;/2TL) 

X Im[xi(g,a;,TL)]Im[xj(g,a;,TL)]. 

This eqnation can be nsed to dehne collective electronic cooling times which are easier to 
compare to the experimental carrier relaxation times extracted from the phenomenological 
hts. Fignre compares the experimental and the theoretical relaxation times over a lattice 
temperatnre range of Tl = 10 — 160 K. From this hgnre we see that, within onr approxima¬ 
tions, theory reprodnces the experimental trends versns both the lattice temperatnre and the 
nnmber of epitaxial graphene layers; more specihcally, both relaxation times increase with 
decreasing lattice temperature and increasing the number of layers. We note that since our 
theoretical relaxation times are obtained by linearizing the interlayer energy transfer rate 
equations, they underestimate interlayer coupling except when Tld — Thd "C Thd- Thus, we 
expect (and observe) the theoretical relaxation times to be longer than those obtained by 
htting the experimental data over a broader temperature range. 

We can identify the microscopic origins of the relaxation times’ dependence on lattice 
temperature using a formula derived in Supplementary Note 2, where we give an analytic 
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formula for the interlayer Coulombic energy transfer rate between a single LD layer and a 
single HD layer in MEG. Here, we show only the result of summing over the rates between all 
HD-LD pairs of layers. Applying this result to the multilayer case, we obtain the following 
net energy loss rate (per area) for N LD layers in MEG: 


^ • Qhd,ld =7(Tld - Ti 


ATz/, 


LD 


L2 


HDJ 


(Zl,- ^^hd) 




T, 


F,LD 


T, 


HD 


(13) 


where 7 = /30){k^/h^Vp) and the density of states in the HD and LD layers are denoted 

by r'HD and i^ld- The sum over index j runs over all HD layers in the MEG sample. We 
hnd that the interlayer Goulombic energy transfer rate exhibits a temperature dependence 
of T^ln(T).^ This result is related to the familiar result^ for electron-electron scattering 
rates in Fermi liquids which are proportional to T2ln(T). These power laws appear, 
because the Fermi distribution limits the initial-hnal state pair energies which can partake 
in scattering to an energy window of width hu ~ k^T about the Fermi energy, which 
qualitatively explains the additional factor of T; the additional factor of T appears because 
scattering events are weighted by the transferred energy in the cooling power case. On the 
other hand, the heat capacity of nearly neutral graphene increases linearly with T (or as if 
T Tp). The hnal result is that the electronic cooling time dehned above becomes shorter 
with increasing lattice temperature, consistent with the experimental trends as shown in 
Figure Precise agreement is not expected outside of the degenerate temperature regime 


(Tl -C Tf,ld ~ 135 K), within which Equation 13 becomes exact 


Discussion 

In conclusion, we have developed a theory of hot-carrier equilibration based on interlayer 
energy transfer via screened Goulomb interactions between electrons in the many top low- 
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density (LD) layers and in the few high-density (HD) layers close to the nnderlying SiC 
snbstrate in mnltilayer epitaxial graphene (MEG). The theory is complicated by the essential 
role of dynamic screening of the electron-electron interactions in all MEG layers throngh 
the temperatnre-dependent charge carrier response. To obtain a transparent theory, we 
have made two well-jnstihed simplihcations. First, we note that screening in the relevant 
temperatnre, wave vector, and freqnency regime is dominated by the hrst few HD layers close 
to the snbstrate, allowing ns to neglect the screening by the top LD layers. Second, we note 
that interlayer energy transfer among the LD layers is mnch stronger than between LD and 
HD layers, and we therefore can describe all LD layers by a common electron temperatnre. 
We have compared the calcnlated cooling dynamics with the relaxation dynamics measnred 
via nltrafast time-resolved THz spectroscopy. The observed experimental dynamics exhibit 
the expected timescales, dependence on lattice temperatnre, and dependence on nnmber of 
epitaxial graphene layers predicted by the theory, providing strong snpport for the proposed 
mechanism, within the approximations necessary in the development of the theory. The 
theoretical approach developed here may be expected to be applicable to many other types 
of layered 2D electron systems. These may inclnde semicondnctor heterostrnctnres as well 
as the wide variety of novel 2D materials nnder active development inclnding transition 
metal dichalcogenides (e.g. M 0 S 2 , MoSe 2 , WS 2 , WSe 2 , etc.),® and other van der Waals 
heterostrnctnres.®^ 
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Figure 1: Interlayer Coulombic energy transfer between two graphene layers a, 

Schematic diagram of interlayer Coulombic energy transfer from a hot LD to a cold HD 
graphene layer, b, Cooling power of interlayer Coulombic energy transfer as a function 
of electron temperature for Fermi level F^f,hd = 300 meV in the HD graphene layer and 
various Fermi levels F^f,ld in the LD graphene layer, c, Ratio of the cooling power of in¬ 
terlayer Coulombic energy transfer to the cooling power of intralayer acoustic phonon 
cooling as a function of electron temperature for i?F,HD = 300 meV and various values 
of i?F,LD showing that interlayer Coulombic energy transfer can dominate intralayer acous¬ 
tic phonon cooling in the LD graphene layer, d, Ratio of the cooling power of interlayer 
Coulombic energy transfer Q®* to the cooling power of intralayer disorder-assisted electron- 
phonon (supercollision) cooling Q®® as a function of electron temperature for F^f,hd = 300 
meV, F^f.ld = 10 meV (typical for C-face MEG on SiC) and various values of the disorder 
mean free path in the LD graphene layer showing that interlayer Coulombic energy transfer 
can dominate intralayer disorder-assisted electron-phonon (supercollision) cooling in the LD 
graphene layer. 
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Figure 2: Ultrafast time-resolved THz spectroscopy on MEG a, Schematic diagram 
of the ultrafast time-resolved THz spectroscopy set-up. (Inset: Schematic diagram of a MEG 
sample with a gradient doping density prohle.) b, Normalized differential THz transmission 
spectra At{u)/t{u) recorded at a pump ffuence of 0.87 /rJ cm“^ and a substrate temperature 
of 40 K for a few different pump-probe delays for a MEG sample with ~ 63 layers. The 
black dashed line indicates the experimental noise level. The THz spectra are remarkably 
dispersionless in the detectable frequency range under all experimental conditions, c-d, 
Linear (c) and logarithmic (d) plots of normalized differential THz transmission at the peak 
of the THz probe pulse At/t as a function of pump-probe delay recorded at a pump ffuence 
of 23.4 yuJ cm“^ for a few different substrate temperatures for a MEG sample with ~ 63 
layers. The THz carrier dynamics evolve from a faster mono-exponential relaxation at room 
temperature to a slower bi-exponential relaxation at cryogenic temperatures. Subfigures (c) 
and (d) share the same legend, e-f, Short and long relaxation times as a function of substrate 
temperature for a few different pump fiuences for a MEG sample with ~ 63 (e) and 35 (f) 
layers. The values are extracted from phenomenological fits to normalized differential THz 
transmission At/t. The long relaxation times increase with the number of epitaxial graphene 
layers, which indicates the presence of interlayer interaction in MEG with HD and LD layers, 
g, Normalized differential THz transmission at the peak of the THz probe pulse At/t as a 
function of pump-probe delay recorded at a pump ffuence of 60.0 /xJ cm“^ for a few different 
substrate temperatures for a MEG sample with ~ 3 layers. The THz carrier dynamics 
follow a fast mono-exponential relaxation at all temperatures, h, Relaxation times as a 
function of substrate temperature for a few different pump fiuences for a MEG sample with 
~ 3 layers. The relaxation times are completely independent of the substrate temperature, 
because there is practically very little or no interlayer energy transfer in MEG with all HD 
layers. The experimental error in all relaxation times is due primarily to long-term drift of 
the optomechanical components and the ultr^st Ti:Sapphire laser system, and is estimated 
not to exceed 
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Figure 3: Ultrafast degenerate IR pump-probe spectroscopy on MEG Normalized 
probe differential transmission AT/T in ultrafast degenerate 1.8 /im IR pump-probe spec¬ 
troscopy as a function of pump-probe delay recorded at a pump fluence of 80 /iJ cm“^ and 
a substrate temperature of 10 K for a MEG sample with ~ 63 layers. The black solid line 
is a guide for the eye. The sign changes in the differential transmission at ~ 1 and ~ 20 ps 
indicate the presence of interlayer energy transfer from the top layers to the first HD layer 
in MEG. (Inset: Schematic diagram of interlayer Goulombic energy transfer from a hot LD 
to a cold HD layer in MEG. The pump selectively injects hot electrons in all layers of MEG 
except the hrst HD layer, in which interband absorption is Pauli blocked.) 
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Figure 4: Interlayer Coulombic energy transfer in MEG Collective cooling dynamics 
of the LD layers in MEG with 30 and 60 layers, calculated using Equation Subfigure (a) 
is for lattice temperature Tl = 10 K and Subfigure (b) is for lattice temperature Tl = 50 
K. Both figures illustrate that the MEG hot-carrier relaxation characteristics observed in 
Figure are naturally described by the interlayer Goulombic energy transfer mechanism. 
Samples with more LD layers cool more slowly because Goulomb coupling to the HD layers 
close to the substrate decreases with increasing layer separation. 
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Figure 5: Comparison between experimental and theoretical relaxation times in 

MEG Subfigure (a) shows experimental hot-carrier relaxation times versus substrate tem¬ 
perature for MEG samples with ~ 63 and ~ 35 layers. The values correspond to the long 
relaxation times from Figure |^-f. Subfigure (b) shows theoretical hot-carrier relaxation 
times versus substrate temperature calculated using the linearized interlayer Coulombic en¬ 
ergy transfer rate Equation 12 The interlayer Coulombic energy transfer mechanism ex¬ 


plains relaxation time trends versus substrate temperature and multilayer system thickness. 
The factor of two discrepancy in magnitude is typical for RPA theories of electron-electron 
scattering amplitudes. The level of agreement found here is similar to that found in compar¬ 
isons of RPA calculations and measurements of the Coulomb-coupled interlayer momentum 
transfer rate (i.e. Coulomb drag resistivity) between neighboring GaAs/AlGaAs quantum 
wells.^ 
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Supplementary Information 


Supplementary Note 1 

Normalized differential THz transmission spectra At{uj)/t{Lu) 

Here, we present additional experimental data to further demonstrate that the normal¬ 
ized differential THz transmission spectra are remarkably dispersionless in the detectable 
frequency range under all experimental conditions. Supplementary Figure [T] and Supple¬ 
mentary Figure show the differential THz transmission spectra normalized to the THz 
transmission without photoexcitation, At{u)/t{u), for variable substrate temperature and 
variable pump fluence, respectively, for a MEG sample with ~ 63 layers. The frequency- 
independent dynamic THz response of the MEG samples justifies recording the normalized 
differential THz transmission only at the peak of the THz probe pulse, At/t, as a function 
of pump-probe delay to map out the relaxation dynamics of the photoexcited carriers. The 
slight fluctuations in a few of the data scans for frequencies >1.7 THz are due to water va¬ 
por absorption arising from very slight fluctuations in the humidity level between the sample 
and the reference scans. We note that spectra obtained from Fourier-transformed time- 
domain measurements are much more susceptible to noise than direct spectrally-resolved 
measurements, because (long-term) fluctuations in the time domain are transferred into un¬ 
certainties in the frequency domain during the Fourier transformation process. On the other 
hand, direct time-resolved THz spectroscopy measurements of the relaxation dynamics are 
less prone to noise from fluctuations, and a very high signal-to-noise ratio (SNR) can be 
achieved through sufficient integration. 
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Supplementary Note 2 

Interlayer energy transfer linearized in AT 

Here, we derive an approximation for the interlayer Coulombic energy transfer rate between 
a lightly doped (LD) layer and a highly doped (HD) layer in MEG systems. We neglect 
electron tunneling between layers. Energy transfer between layers can nevertheless occur 
when an electron in one layer scatters off an electron in a remote layer. To hrst order in 
the temperature difference ST = — Thd, we find that the energy transfer rate between 

thermal electron distributions in two Coulomb-coupled layers is given by: 



4Tl sinh^(ha;/2THD) 

X lm[xi{q, cu, THD)]Im[xj(g, u, Thd)]- 


Equation can be derived by summing over all interlayer electron collision processes and 
combining a Fermi golden-rule expression for the transition rates with a random phase ap¬ 
proximation expression for the electron-electron scattering amplitudes. We focus on the 
low temperature limit (i.e. T ^ Tf^ld,Tf^hd), where several physical approximations can 
be made. In this case, it is clear that the a; —>■ 0 limit dominates the integrand of Equa¬ 
tion We note that in this low frequency limit, Im[xj(g,ci; —)■ 0)] = —h>iOj/{v-pq)^ where 
z/j = 2TF,j/(7rh^n|) is the density of states at the Fermi energy. This linear dependence on 
frequency is also found in parabolic band 2DEG’s, and describes how the particle-hole exci¬ 
tation spectrum vanishes with decreasing excitation energy fujj^ Additionally, in the degen¬ 
erate regime interlayer particle scattering allows a maximum change in electron wavevector 
of Agmax = 2 A;f, a fact reflected in: 

Im[xi(g,a;-)■ 0)] = 0, g > 2 /cf- (2) 
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The interlayer Coulomb interaction, proportional to naturally places the limit g < l/c?, 
where d is the interlayer separation. However, /cf.hd ^ ^/d and we can approximately 
neglect interlayer separation (i.e. d —?• 0) in Equation Then, the MEG dielectric function 
reduces to a Thomas-Fermi-like form: 


MEG _ 

4EG(,) = 1 + ,meg ^ ^ ^ ,3) 

where uj is the density of states in the j’th HD layer and k is the background dielectric 
function of a thin him on SiC {k ~ (10 + l)/2 = 5.5). Because of the large difference in 
carrier density between layers near and far from the substrate, allowed q values are much 
smaller than and the screened interlayer interaction reduces to (27re^)/(«:g^^'^). 

The remaining integrals in Equation [T] contain a logarithmic divergence at zero wavevec- 
tor. This is removed using a cutoff of uj/v-p which rehects the zero value of Im[y;(g,ci;)] above 
the intraband particle-hole continuum. 

Finally, we hnd that the linearized interlayer energy transfer rate per area between a pair 
of graphene layers in a MEG system (where one is HD and the other LD) is: 


1 

Z2 


• Q 


el 

HD,LD 


=7(^ld 


Tub) 


T" In 







( 4 ) 


where 7 = ( 87 r^/ 30 )(fc|/h^np). The density of states in the HD layer (LD layer) is denoted by 
^HD (j^ld)- The sum over index j runs over all HD layers in the MEG system, accounting for 
the approximation that the dominant screening effect originates from the majority fraction 
of carriers in the HD layers (see the main text). Tppo = T^f,ld/^b is the Fermi temperature 
of the LD layer. Equation becomes exact in the degenerate limit T <C Tf,ld,Tf,hd- 
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Supplementary Note 3 

Interlayer energy transfer to leading order in T 

Here, we derive in detail the leading order in temperature formula for the interlayer Coulom- 
bic energy transfer rate between one lightly doped (LD) layer of graphene and one highly 
doped (HD) layer of graphene. This follows similar steps as in the previous section, but 
without the assumption of inhnitesimal temperature separation between the LD and HD 
layers. As described in the main text, the interlayer energy transfer rate per area between 
an LD and an HD layer is: 


Qel 

L2 


h 


47^3 


udu dg I^LD hoT [’^b(^ld) —''^b(^hd). 


X Im[xLD(g, UJ, TLD)]Im[xHD(g, UJ, Thd)], 


(5) 


where nB(T) is a Bose distribution function and Uldhd = is the 

screened interaction between electrons in opposing layers within the random phase approxi¬ 
mation (RPA). Given that the separation of the two layers is a distance d, the RPA dielectric 
function is given by: 


—(1 — VqXLD{q,<^,TBD)){^ — VqXHD{q, 

— VgC ^^'^XLD(g,l^,dLD)XHD(g, i^jThd), 


( 6 ) 


where Vq = 2ne^/nq and Xj(g, a;,T) is the temperature dependent non-interacting density- 
response function of the i’th layer of graphene. We next assume the temperature of the 
electrons in the HD layer is pinned to the lattice temperature, and approximate this as zero 
relative to the high temperature in the LD layer, i.e. Thd = Tl —t 0. Relabeling Tld as T 
we have: 


Qel 

L2 


h 



dq 


VqC 


—qd 




'nB(T) 


X Im[xLD(g,ca,T)]Im[xHD(g,ca)L 


(7) 
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In the degenerate limit, k-^T « -Ef,ld, the Bose distribution function limits the important 
frequencies to approximately u < -Ef.ld/ h, and reveals that it is the leading order in u which 
will yield the leading order in T. With this motivation we make use of the low frequency 
limit, Im[xj(g,a; —)■ 0)] = —i'iUj/{vYq)i where z/j = 2i?F,i/(7rh^nf) is the density of states at 
the Fermi energy. Similarly, the dielectric function can be reduced to: 


,RPA/ 


TF 

9ld 


e"*‘'(g,a;,T) —)■ (1 H-)(1 + 

q 


TF 

feO' 


— e 


TF TF 
-2(?d?LD?HD 


( 8 ) 


where the Thomas-Fermi wavevector of the Fth layer is dehned as qj^ = qVqUi. Making use 
of the fact that in our MEG samples /cnpod << 1 and 5 'hd/^f,ld >> 1, we simplify the 
interlayer energy transfer rate per area to: 


Qcl 

L2 


-^F.LD^LD 


VtdVt / QdQ 




(9) 


where we have introduced the dimensionless wavevector, Q = g//cF,LD 5 the dimensionless fre¬ 
quency, G = tvjj/E y;lDi and the dimensionless temperature, t = k^T/E- p;ld- The wavevector 
integration here diverges logarithmically. To remedy this, we identify that Im[xLD(Q' 5 0)] 
vanishes for uj > v^q, which cuts off the wavevector integration for Q < Q. Carrying out the 
remaining integrals, we obtain the leading order in temperature energy loss rate per area of 
the LD electrons: 


L2 


( 10 ) 
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Supplementary Note 4 

Interlayer energy transfer with no LD-LD layer coupling 

Here, we illustrate the distance dependence of the interlayer Coulombic energy transfer rate 
by calculating the temperature dynamics of an individual LD layer coupled to the HD layers 
in MEG systems, when we neglect energy transfer between pairs of LD layers. Although 
Coulomb coupling between LD layers is very strong,® by temporarily forcing ld' 0 
we can gain insight into the distance dependence of the interlayer energy transfer rate. 
The temperature dynamics of each individual LD layer coupled to the HD layers are then 
independent and obey: 



( 11 ) 


If we also approximate the heat capacity in the LD layers by the neutral graphene formula 
Cld = 18C(3)^ld/(^'^f)) obtain the results shown in Supplementary Figure where we 
have used this formula to evaluate Tu:,{t) for several different values of (iHD,LD- All of these 


curves are calculated for the typical carrier density of the LD layers in MEG, measured in 
experiment, of uld = 10^°cm“^. The slowest (dHD,LD = 60 layers) and fastest ((iHD,LD = 1 
layer) temperature relaxation curves provide upper and lower bounds, respectively, on the 
true relaxation time of all LD layers when interlayer energy transfer between pairs of LD 
layers is no longer neglected, i.e. Qld ld' 7^ 0- 
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Supplementary Note 5 

Acoustic phonon pinning of Thd to Tl 

Here, we present some simple calculations in support of our assumption that acoustic phonon 
cooling is capable of pinning the electronic temperature in the HD layers at the lattice tem¬ 
perature while these layers act as a heat sink for energy dissipation from the remaining hot 
LD layers. Heuristically, we also note that previous ultrafast optical spectroscopy experi¬ 
ment^ have observed thermal relaxation times in the HD layers of MEG of ~ 10 ps, much 
faster than the relaxation times in the LD layers of ~ 100 — 500 ps that we report here. 

Acoustic phonon cooling has previously been investigated in the context of hot-carrier 
cooling in disorder-free monolayer graphene.^El In these single layer systems, as a result of 
the relatively large optical phonon energy in graphene {hwop ~ 200 meV^, acoustic phonons 
serve as the primary intrinsic cooling mechanism over a large temperature window extending 
up towards T ~ 250 K.^l To estimate the ability of acoustic phonon cooling to take away 
the electronic energy that is Coulomb-transfered from LD to HD layers, we use Equation 14 
in Bistritzer and MacDonald® to calculate the total energy transfer rate to the the lattice. 
Using the carrier density prohle of the four most highly doped layers measured in experiment 
{Ep^i _4 = 360,218,140, and 93 meV), we hnd that = 9.09(Thd — Tp) W cm“^ can be 
transferred to the lattice via carrier-phonon scattering in the HD layers. We can calculate 
the ratio of to the interlayer energy transfer rate from the N LD layers to these 
four HD layers using Equation 13 in the main text. Supplementary Figure suggests that 
for both iV = 30 and iV = 60 layer MEG systems, the acoustic phonon cooling power is 
sufficient to keep the HD layers pinned to the lattice temperature while absorbing energy 
from the distant LD layers. 
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Supplementary Note 6 

Disorder-assisted electron-phonon (supercollision) cooling 

Here, we present an application of the recently proposed disorder-assisted electron-phonon 
(supercollision) cooling mechanisml^ to MEG. We investigate the qualitative and quantitative 
differences between disorder-assisted electron-phonon cooling in HD and LD graphene and 
we clearly demonstrate that this cooling mechanism cannot alone explain electronic cooling 
in high quality MEG. 

The electron temperature dynamics of a single graphene layer in the framework of the 
disorder-assisted electron-phonon cooling mechanism is given by: 

dtT = Q^yc, (12) 

where = dt£ is the electronic cooling rate due to supercollisions, C = drS is the electronic 
heat capacity and S is the electronic energy density. The electronic cooling rate due to 
supercollisions in the degenerate limit when k-sT Ep is given by:® 

= -A{T^ -Tl), (13) 

with a rate coefficient A = 'd.Q2g^v'^{E-p)k^/{hk-pl), where g = D/\/2pv'^ is the electron- 
phonon coupling constant and i^(Ef) = EF/{27rh‘^Vp) is the density of states at the Fermi 
level per one spin and one valley flavor. The rest of the parameters are the Fermi velocity 
Up, the sound velocity Vs, the deformation potential D, the mass density p and the disorder 
mean free path 1. The electronic cooling rate in the non-degenerate limit when k^T 3> Ep 
is given by:® 

= -B{T^ -Tl), (14) 

with a rate coefficient B = {{4kp)/{Ep))A. Because both the supercollision cooling rate and 
the electronic heat capacity have different functional dependence on the electron temperature 
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at high and at low doping densities, we need to consider the two cases separately. 


For HD graphene, the heat capacity is given by C = aT = {{27iEFk'^)/{3h?Vp))T. By 
snbstitnting this expression and Eqnation 13 in Eqnation 12, we obtain:® 


AT^ 

rri 

a 1 


(15) 


The electronic cooling timescale is governed by the rate coefficient A/a <x. E-p/k-pl oc 1//. In 


the high and low electron temperatnre limits, Eqnation 15 rednces to: 


d,T = --T^ for T>Tl, 
a 


3A 

a,T =- Ti,(T - Tl) for T « Tl, 

a 


(16) 

(17) 


with solntions given by: 


T{t) = 

1 + int 

for 

T>Tl, 

(18) 

T{t)=n + {To-n)exp( 

—nt 

a / 

for T oi Tl. 

(19) 


For LD graphene, the heat capacity is approximated by C = = ((18C(3)fc|)/(7r/i^n|))T^. 

By snbstitnting this expression and Eqnation in Eqnation we obtain: 


BT^ -Tl 


(3 T2 


( 20 ) 


The electronic cooling timescale is governed by the rate coefficient B/(3 oc l/kpl oc 1/Epl. 
In the high and low electron temperatnre limits, Eqnation rednces to: 


dtT 



for T > Tl, 


( 21 ) 
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( 22 ) 


aT = 


fTliT-T,) 


with solutions given by: 


for T ^ Tl, 


T(t) = 


Tn 


l + 2fT2t 


/or T > Tl, 


T(t) = Tl + (To - Tl) exp ( -y T/t 


5T. 


for 


T ^ Tl. 


(23) 


(24) 


We observe that the different functional form of the cooling rate and the heat capacity 
at high and at low doping densities results in qualitatively different electron temperature 
relaxation dynamics. In the low electron temperature limit, in particular, the temperature 
dynamics follows an exponential form with a lattice temperature dependent characteristic 
time thd = a)Tif)~^ oc //Tl for HD graphene, and with a lattice temperature dependent 
characteristic time tld = ((5T//3)Tl)“^ oc E-pl/Tf for LD graphene. 

To obtain quantitative estimates of the full electron temperature dynamics predicted by 


the disorder-assisted electron-phonon cooling mechanism, we solve Equation 15 and Equa¬ 
tion 2^ numerically. In all calculations, we use up = 1 x 10® m s“^, Vg = 2.1 x 10^ m s“^, 
p = 7.6 X 10“^ kg m~^ and T = 20 meV. Supplementary Figure]^ and Supplementary Fig¬ 
ure 1^ show the calculated electron temperature dynamics for HD graphene with E-p = 100 
meV and for LD graphene with Ep = 10 meV, respectively, at Tl = 10 K for variable disor¬ 
der mean free path /. The large uncertainty in the value of the disorder length scale leads 
to a very large spread in the calculated electron temperature dynamics. 

The experimental measurement of the precise value of the disorder mean free path be¬ 
tween supercollisions is a challenging task and reliable estimates are lacking in the literature. 
One reliable way to characterize the degree of disorder which can potentially be related to 
supercollision cooling is directly from the width of the Dirac cone near the Dirac point in 
the graphene band structure as measured in high-resolution angle-resolved photoemission 
spectroscopy (ARPES). We use such recent ARPES measurements to estimate the disorder 
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mean free path value appropriate for our MEG samples.^ In Ref. [TU] (Figure SI), correlation 
lengths of ~ 1 — 3 nm in exfoliated graphene and correlation lengths exceeding ~ 50 nm 
(limited only by the instrument resolution, but are expected to be even longer) in C-face 
MEG have been reported. Supplementary Figure shows the calculated electron temper¬ 
ature dynamics for HD graphene with Ep = 100 meV and disorder mean free path I = 5 
nm for variable Tp. We note that these calculations are consistent with recent experimen¬ 
tal results based on photocurrent measurements on HD chemical-vapor-deposited (GVD) 
graphene.® Supplementary Figure]^ shows the calculated electron temperature dynamics for 
LD graphene with Ep = 10 meV and disorder mean free path I = 50 nm for variable Tp. 

We also present the electronic cooling times in the low electron temperature limit pre¬ 
dicted by the disorder-assisted electron-phonon cooling mechanism that are more straightfor¬ 


show the calculated electronic cooling times thd for HD graphene with Ep = 100 meV and 
tld for LD graphene with Ep = 10 meV, respectively, as a function of Tp for variable dis¬ 
order mean free path 1. We observe that the predictions of the supercollision cooling model 
applied to the HD layers of MEG are inconsistent both in magnitude and lattice temper¬ 
ature dependence with the ultrafast time-resolved THz spectroscopy experiments. We also 
observe that the model predictions for the LD layers of MEG can roughly capture the order 
of magnitude and the lattice temperature dependence for some disorder mean free path, but 
not the layer number dependence. Because the quality of our MEG samples is expected to 
be even higher than we have conservatively assumed in these calculations, we conclude that 
the disorder-assisted electron-phonon (supercollision) cooling can provide a parallel cooling 
channel, but it is not dominant in MEG. 


ward to compare to the experiment. Supplementary Figured and Supplementary Figure 10 
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Supplementary Figure 1; Ultrafast time-resolved THz spectroscopy on MEG. Nor¬ 
malized differential THz transmission spectra At{u)/t{u) recorded at a pump fluence of 0.87 
fiJ cm“^ and a pump-probe delay of 1 ps for a few different substrate temperatures for a 
MEG sample with ~ 63 layers. The black dashed line indicates the experimental noise level. 



Frequency (THz) 

Supplementary Figure 2: Ultrafast time-resolved THz spectroscopy on MEG. Nor¬ 
malized differential THz transmission spectra At{uj)/t{uj) recorded at a substrate tempera¬ 
ture of 280 K and a pump-probe delay of 1 ps for a few different pump fluences for a MEG 
sample with ~ 63 layers. The black dashed line indicates the experimental noise level. 
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Supplementary Figure 3: Interlayer energy transfer with no LD-LD layer coupling. 

Temperature dynamics when interlayer energy transfer between LD layers (uld = 10^°cm“^) 
is ignored. The LD layer temperature Tld(^) resulting from cooling via interlayer Coulombic 
energy transfer to the HD layers near the substrate (uhd ^ lO^^cm”^) at a constant lattice 
temperature Tl. The distance between the particular LD layer and the HD layers is varied, 
duD.LD = 60a, 50a, 40a, 30a, 20a, 10a, a (top to bottom) where a = 3.4 Angstroms. Subfigure 
(a) shows results for lattice temperature Tl = 10 K and Subfigure (b) for Tl = 50 K. 
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Supplementary Figure 4: Acoustic phonon pinning of Thd to Tl. Ratio of the acoustic 
phonon cooling power in the HD layers of MEG to the interlayer Coulombic energy trans¬ 
fer rate from the LD to the HD layers of MEG as a function of the lattice temperature 
Tl. 
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Supplementary Figure 5: Disorder-assisted electron-phonon (snpercollision) cooling 
in HD graphene. Electron temperature dynamics T{t) — Ti^ predicted by the disorder- 
assisted electron-phonon cooling mechanism for HD graphene with E-p = 100 meV at Tp = 10 
K for variable disorder mean free path 1. 



Supplementary Figure 6: Disorder-assisted electron-phonon (snpercollision) cooling 
in LD graphene. Electron temperature dynamics T(t) — Tp predicted by the disorder- 
assisted electron-phonon cooling mechanism for LD graphene with Ep = 10 meV at Tp = 10 
K for variable disorder mean free path 1. 
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Supplementary Figure 7; Disorder-assisted electron-phonon (snpercollision) cooling 
in HD graphene. Electron temperature dynamics T{t) — Ti^ predicted by the disorder- 
assisted electron-phonon cooling mechanism for HD graphene with E-p = 100 meV and 
disorder mean free path I = 5 nm for variable Tp. 



Supplementary Figure 8: Disorder-assisted electron-phonon (snpercollision) cooling 
in LD graphene. Electron temperature dynamics T(t) — Tp predicted by the disorder- 
assisted electron-phonon cooling mechanism for LD graphene with Ep = 10 meV and disorder 
mean free path I = 50 nm for variable Tp. 
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Supplementary Figure 9; Disorder-assisted electron-phonon (snpercollision) cooling 
in HD graphene. Electronic cooling time thd in the low electron temperature limit pre¬ 
dicted by the disorder-assisted electron-phonon cooling mechanism for HD graphene with 
Ey = 100 meV as a function of Tl for variable disorder mean free path 1. 



Supplementary Figure 10: Disorder-assisted electron-phonon (snpercollision) cool¬ 
ing in LD graphene. Electronic cooling time tld in the low electron temperature limit 
predicted by the disorder-assisted electron-phonon cooling mechanism for LD graphene with 
Ep = 10 meV as a function of Tl for variable disorder mean free path 1. 
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